Example 2: Currency Substitution with capital accumulation. See paper for details. (Corresponds to Example 3 of the longer version of the paper.
Manoj Atolia: Septermber 18, 2007.
Works with both Mathematica 6.0 and 5.2.

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[FindMinimum::lstol];Off[FindMinimum::reged];Off[InterpolatingFunction::dmval];Off[Precision::precsm];

(* The precision of preliminary calculations. *)
{precision = 100, $MinPrecision=digits=precision};

Initialization

(* Set the directory where Mathematica will store intermediate results. *)
SetDirectory["c:"];

(* ALWAYS ENTER NUMERICAL VALUES AS FRACTIONS AND NOT AS DECIMALS TO ENSURE REQUISITE PRECISION IN CALCULATIONS *)   
{τ=(1/4),ρ=(1/10),σ=2,δ=(1/20),β=(1/2),Ω=2}//N                                     (* Parameters of the model *)

{(Subscript[bbyY, o])=(25/100),(Subscript[mbyF, o])=(1/2),(Subscript[mbyY, o])=(15/200),(Subscript[π, o])=(1/2),(Subscript[θ, o])=(1/2)}//N           (* Data for calibration *)

{(Subscript[z, o])=1,(Subscript[k, o])=1}//N                                                                                                    (* Normalization *)

{(Subscript[z, f])=(6/5)(Subscript[z, o]),ω=(1/5)}//N                                                                                            (* Shock - Exogenous rise in productivity *)

Calibration and Steady State

(* Define utility function, production function etc. *)
U[c_]:=(1/(1-(1/τ)))(c^(1-(1/τ)));(Subscript[U, 1])[c_]:=Evaluate[D[U[c],c]];ϕ[m_,F_]:=((Subscript[g, 0])/(1-1/τ))(((((((Subscript[g, 1])(m^((σ-1)/σ))+(1-(Subscript[g, 1]))(F^((σ-1)/σ))))^(σ/(σ-1)))))^(1-1/τ));(Subscript[ϕ, 1])[m_,F_]:=Evaluate[D[ϕ[m,F],m]];(Subscript[ϕ, 2])[m_,F_]:=Evaluate[D[ϕ[m,F],F]];
Y[k_,z_]:=z((((Subscript[a, 0])+(Subscript[a, 1])(k^((β-1)/β))))^(β/(β-1)));(Subscript[Y, 1])[k_,z_]=Evaluate[D[Y[k,z],k]];
(* Avoid using symbols for subscripts for representing derivatives. Thus, (Subscript[Y, 1])[k,z] is (Subscript[Y, K])[k,z] of the paper.*)
ψ[x_]:=(n/2)(((x-δ))^2);n=(1/(Ω δ));
(* Since we are using symbols o and f as subscripts for initial and final steady states avoid using o and f as names of variables or function. *)

(* Calibrate and compute the original steady state *)
{(Subscript[Y, o])=Y[(Subscript[k, o]),(Subscript[z, o])];,(Subscript[c, o])=(Subscript[Y, o])-δ (Subscript[k, o]),(Subscript[m, o])= (Subscript[mbyY, o]) (Subscript[Y, o]),(Subscript[F, o])=((Subscript[m, o])/(Subscript[mbyF, o])),(Subscript[b, o])=(Subscript[bbyY, o]) (Subscript[Y, o]),L=(Subscript[π, o])(Subscript[m, o])-ρ (Subscript[b, o])};
{(Subscript[g, 0])=.,(Subscript[g, 1])=.};eq1:= (Subscript[ϕ, 1])[(Subscript[m, o]),(Subscript[F, o])]== (Subscript[U, 1])[(Subscript[c, o])](ρ+(Subscript[π, o]));eq2:=(Subscript[ϕ, 2])[(Subscript[m, o]),(Subscript[F, o])]== (Subscript[U, 1])[(Subscript[c, o])]ρ ;{(Subscript[a, 0]),(Subscript[a, 1])}=.;{(Subscript[a, 0]),(Subscript[a, 1])}={(Subscript[a, 0]),(Subscript[a, 1])}/.FindRoot[{(Subscript[Y, 1])[(Subscript[k, o]),(Subscript[z, o])] == ρ+δ,(Subscript[θ, o])==(((ρ+δ)(Subscript[k, o]))/(Subscript[Y, o]))},{(Subscript[a, 0]),1},{(Subscript[a, 1]),1},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];
{(Subscript[g, 0]),(Subscript[g, 1])}={(Subscript[g, 0]),(Subscript[g, 1])} /. FindRoot[{eq1,eq2},{(Subscript[g, 0]),(10^(-6))},{(Subscript[g, 1]),1},MaxIterations->1000,AccuracyGoal->$MinPrecision,WorkingPrecision->digits];

(* Compute the new steady state *)
(Subscript[k, f])=.;(Subscript[k, f])=(Subscript[k, f])/.FindRoot[{(Subscript[Y, 1])[(Subscript[k, f]),(Subscript[z, f])]==ρ+δ},{(Subscript[k, f]),(Subscript[k, o])},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits];(Subscript[Y, f])=Y[(Subscript[k, f]),(Subscript[z, f])];(Subscript[c, f])=(Subscript[Y, f])-δ (Subscript[k, f]);{(Subscript[m, f]),(Subscript[F, f]),(Subscript[π, f])}=.;{(Subscript[m, f]),(Subscript[F, f]),(Subscript[π, f])}={(Subscript[m, f]),(Subscript[F, f]),(Subscript[π, f])}/.FindRoot[{(((Subscript[ϕ, 2])[(Subscript[m, f]),(Subscript[F, f])])/((Subscript[U, 1])[(Subscript[c, f])]))==ρ,(((Subscript[ϕ, 1])[(Subscript[m, f]),(Subscript[F, f])])/((Subscript[U, 1])[(Subscript[c, f])]))==ρ+(Subscript[π, f]),(Subscript[π, f])(Subscript[m, f])==L+ρ (Subscript[b, o]) },{(Subscript[m, f]),(Subscript[m, o])},{(Subscript[π, f]),(Subscript[π, o])},{(Subscript[F, f]),(Subscript[F, o])},AccuracyGoal->$MinPrecision,WorkingPrecision->digits];
{{(Subscript[m, o]),(Subscript[F, o]),(Subscript[k, o]),(Subscript[c, o]),(Subscript[Y, o]),(Subscript[b, o]),L,(Subscript[z, o]),(Subscript[π, o]),(Subscript[c, o])/(Subscript[m, o]),(Subscript[m, o])/(Subscript[Y, o]),(Subscript[m, o])/(Subscript[F, o])},{(Subscript[m, f]),(Subscript[F, f]),(Subscript[k, f]),(Subscript[c, f]),(Subscript[Y, f]),(Subscript[b, o]),L,(Subscript[z, f]),(Subscript[π, f]),(Subscript[c, f])/(Subscript[m, f]),(Subscript[m, f])/(Subscript[Y, f]),(Subscript[m, f])/(Subscript[F, f])}}//N//MatrixForm

Global Nonlinear Saddlepath

(* Specify the nonlinear system with state variables at the end *)
{  mdot = L+r (Subscript[b, o])-pi m;
cdot = τ c (r-ρ);
Fdot = Y[k,z]-c-i-ψ[i/k]k;
kdot=i-δ k;
zdot =  ω((Subscript[z, f])-z);
eq0:= ((ψ''[i/k])/(1+ψ'[i/k]))(i/k)((idot/i)-(kdot/k))==r+δ-(((Subscript[Y, 1])[k,z]+ψ'[i/k]i/k-ψ[i/k])/(1+ψ'[i/k]));idot=.;
idot = First[idot/.Solve[eq0,{idot}]];
system={mdot,cdot,idot,Fdot,kdot,zdot};};

(* Solve for the endogenous variables as functions of the variables in the system. For consistency maintain same order of variables in arguments for functions of endogenous variables as in the system *)
{  r[m_,c_,F_]:= (((Subscript[ϕ, 2])[m,F])/((Subscript[U, 1])[c]));
pi[m_,c_,F_]:=(((Subscript[ϕ, 1])[m,F]-(Subscript[ϕ, 2])[m,F])/((Subscript[U, 1])[c]));};

(* INITIALIZATION CELL - Initialization for finding the required trajectory. *)
{  T = 20;                        (* solution time (years). BEWARE!!! Setting T larger than necessary can substantially increase computational time.*); 
tol = (10^(-4));                (* desired accurcy of the solutions of the variables in the system *);
precision = 16;      (* starting value of precision of calculations*) ;
delta = (10^(-1))(((((((Subscript[F, f])-(Subscript[F, o])))^2)+((((Subscript[k, f])-(Subscript[k, o])))^2)+((((Subscript[z, f])-(Subscript[z, o])))^2)))^(1/2));     (* size of δ-simplex  or distance of simplex  from the final steady state *)};

?solve

(* MAIN CELL - Solving for the global nonlinear saddlepath. *)
Timing[{gtmin,gminsol}=solve[system,{t},{m,c,i,F,k,z},{(Subscript[m, o]),(Subscript[c, o]),δ (Subscript[k, o]),(Subscript[F, o]),(Subscript[k, o]),(Subscript[z, o])},{(Subscript[m, f]),(Subscript[c, f]),δ (Subscript[k, f]),(Subscript[F, f]),(Subscript[k, f]),(Subscript[z, f])},{{r,{m,c,F}},{pi,{m,c,F}}},{T,tol,precision,delta}]]

(* Plot the final solutions. *)
tplot = T;
mPlot=Plot[{If[t<0,(Subscript[m, o]),m[gtmin-t]]/.gminsol,(Subscript[m, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",m},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
cPlot=Plot[{If[t<0,(Subscript[c, o]),c[gtmin-t]]/.gminsol,(Subscript[c, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",c},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
iPlot=Plot[{If[t<0,δ (Subscript[k, o]),i[gtmin-t]]/.gminsol,δ (Subscript[k, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",i},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
FPlot=Plot[{If[t<0,(Subscript[F, o]),F[gtmin-t]]/.gminsol,(Subscript[F, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",F},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
kPlot=Plot[{If[t<0,(Subscript[k, o]),k[gtmin-t]]/.gminsol,(Subscript[k, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",k},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
zPlot=Plot[{If[t<0,(Subscript[z, o]),z[gtmin-t]]/.gminsol,(Subscript[z, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",z},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
YPlot=Plot[{If[t<0,(Subscript[Y, o]),Y[k[gtmin-t],z[gtmin-t]]]/.gminsol,(Subscript[Y, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",Y},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
rPlot=Plot[{If[t<0,ρ,r[m[gtmin-t],c[gtmin-t],F[gtmin-t]]]/.gminsol,ρ},{t,-tplot/10,tplot},PlotRange->All,AxesOrigin->{-tplot/10,ρ},AxesLabel->{"",r},DisplayFunction->Identity];
πPlot=Plot[{If[t<0,(Subscript[π, o]),pi[m[gtmin-t],c[gtmin-t],F[gtmin-t]]]/.gminsol,(Subscript[π, f])},{t,-tplot/10,tplot},PlotRange->All,AxesOrigin->{-tplot/10,(Subscript[π, f])},AxesLabel->{"",π},DisplayFunction->Identity];
If[$VersionNumber==6.0,
Show[GraphicsGrid[{{mPlot,cPlot,iPlot},{FPlot,kPlot,YPlot},{rPlot,πPlot,zPlot}}]],
Show[GraphicsArray[{{mPlot,cPlot,iPlot},{FPlot,kPlot,YPlot},{rPlot,πPlot,zPlot}}]]]

Labels={"m","c","i","F","k","z","Y","r","pi"};InitialSteadyState={(Subscript[m, o]),(Subscript[c, o]),δ (Subscript[k, o]),(Subscript[F, o]),(Subscript[k, o]),(Subscript[z, o]),(Subscript[Y, o]),ρ,(Subscript[π, o])}//N;FinalSteadyState={(Subscript[m, f]),(Subscript[c, f]),δ (Subscript[k, f]),(Subscript[F, f]),(Subscript[k, f]),(Subscript[z, f]),(Subscript[Y, f]),ρ,(Subscript[π, f])}//N;
(Dat=Table[{tt,m[t],c[t],i[t],F[t],k[t],z[t],Y[k[t],z[t]],r[m[t],c[t],F[t]],pi[m[t],c[t],F[t]]}/.{t-> gtmin-t}/.tt->t/.gminsol,{t,0,T,1/5}])//MatrixForm;
(Dat=Join[{Join[{"Initial Steady State"},InitialSteadyState],Join[{"Final Steady State"},FinalSteadyState],Join[{Time},Labels]},Flatten[Dat,1]])//N//MatrixForm

The solve 3D Package

(* Interface *)
BeginPackage["solve`"];

EndPackage[];

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[FindMinimum::lstol];Off[FindMinimum::reged];Off[InterpolatingFunction::dmval];Off[Precision::precsm];

(* Usage messages for the exported functions and the context itself *)
solve::usage = "solve[sys, {t}, vars, varso, varsf, varsendo, init] solves a system of first-order differential equations. \n
sys   = {xdot, ydot,...} lists the right hand side of the system of first-order differential equations (with state variables at end).\n
vars  = {x, y, ...} is the list of variables in the system in the same order as the differential equations in sys.\n
varso = {(Subscript[x, o]), (Subscript[y, o]), ...} is the list of original steady values of variables in the system in the same order as vars.\n
varsf = {(Subscript[x, f]), (Subscript[y, f]), ...} is the list of final steady values of variables in the system in the same order as vars.\n
varsendo = {{w,{x,y, ...}}, ...} is the list of the endogenous variables {w, ...} where w is function of {x,y, ...} where {x,y,.. } are in the order in which they appear as arguments of w. If no endogenous variables in the system, enter {{,{}}}.\n
init   = {tol, precision, epsilon} sets the required accuracy to tol, starting precision to precision, and radius of initial ϵ-ring to epsilon." ;

(* Error messages for the exported objects *)
{  solve::sysDESys = "You called `1` with first argument `2` (the system of differential equations) that has variables other than those defined in `3` and `4`.";

solve::sysLength = "You called `1` with lengths of `2`, `3`, `4` and `5` not equal.";

solve::vars = "You called `1` with third argument `2` that is not a list of variables.";

solve::origVals = "You called `1` with fourth argument `2` that is not a list of numbers.";

solve::finVals = "You called `1` with fifth argument `2` that is not a list of numbers.";

solve::endoVars = "You called `1` with sixth argument `2` that is not a list of lists of the form {w, {x,y, ...}} where w is an endogenous variable and {x,y, ...} are the variables in the system as defined in `3`.";

solve::initVals = "You called `1` with last argument `2` that is not a list of numbers.";

solve::sysError = "You called `1` with a system that does not have three negative eigenvalues when linearized."};

(* Implementation *)
Begin["solve`Private`"];

(* Definition of the exported functions *)
solve[asys_List,{at_},avars_List,avarso_List,avarsf_List,avarsendo_List,ainit_List]:=Block[{sys,t,vars,varso,varsf,varsendo,init,T,tol,precision,delta,Tmax,firstsim,firstrun,precisionini,system,box,A,sysL,dir,localmin,pmin,epsilon,digits,ntrajini,gmindis,gtmin,gminsol,initsol,ini,dis,tsol,tdis,doneflag},

{sys,{t},vars,varso,varsf,varsendo,init}={asys,{at},avars,avarso,avarsf,avarsendo,ainit};
{T,tol,precision,delta}=init;Tmax=4000;firstsim = 01;firstrun = 01;precisionini=precision;

(* Check input *)
checkInput[];

(* Make the system from input data and create box over which to solve the system *)
system=makeSys[];
box:=Apply[Or,Join[MapThread[First[Map[(#1),{t}]]<((Min[#2,#3])/100)&,{vars,varso,varsf}],MapThread[First[Map[(#1),{t}]]>5 Max[#2,#3]&,{vars,varso,varsf}]]];

(* Solve the linearized system *)
A=linearize[];
sysL=Thread[Map[ToExpression[StringJoin[ToString[#],"'","[",ToString[t],"]"]]&,vars]==-A.(Map[ToExpression[StringJoin[ToString[#],"[",ToString[t],"]"]]&,vars]-varsf),List];
dir = Sign[Take[vars,-3]-Take[varsf,-3]/.Thread[vars->sysLsol[t],List]/.t->1000//N];

(* SOLVE FOR THE GLOBAL NONLINEAR SADDLEPATH *);
While[True,
localmin=0;

Label[begin];
$MinPrecision=digits=precision;ntrajini=20;gmindis=tol;

(* Finding the center and the radius of the initial ring for first run of search. *)
If[firstsim==1 && firstrun ==0,{pmin,epsilon}=(<<solution1)[[{3,4}]]];If[firstsim==0 && firstrun ==0, {pmin,epsilon}=(<<solution2)[[{3,4}]]]; 

(* Find the intersection of the linear solution with relevant face of the δ-diamond/simplex to use as the center. *)
If[firstrun ==1,
	epsilon=delta/20;
	
	initsol=FindRoot[dir.Take[sysLsol[t]-varsf,-3]==delta,{t,10},WorkingPrecision->digits];
	pmin=Take[sysLsol[t],-3]/.initsol;Print[{dir,initsol,pmin}];
	
	If[initsol[[1,2]] < T+10  && firstsim ==1 ,
		delta=delta/2;Print["The δ-diamond/simplex is too large. Decreasing 'delta'."];Goto[begin]];
	If[initsol[[1,2]] > T+20  && firstsim==1 ,
		delta = 3delta;Print["The δ-diamond/simplex is unnecessarily tight. Increasing 'delta' by a factor of 3."];Goto[begin]]; ];

(* Conduct an initial check for sufficiency of precision in calculations using the linearized system. *)
If[firstsim==1 && firstrun ==1,
	ini:=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==Apply[Sequence,{stabcond[pmin],pmin},2],List];firstit=1;

	{dis,tdis,tsol}=compdisini[sysL,ini];
	
	Print["LINEARIZED SYSTEM: ","    Minimum Distance from Initial Steady State = ",dis//N,
		"    Time of Closest Approach = ",tdis//N,"    Solution Time = ",tsol//N];
	If[dis <1tol,
		Print["Preliminary check indicates you have enough precision for calculations."],
		precision=precision+4;
		Print["You may not have enough precision for calculations. Increasing 'precision' by 4."];Goto[begin]];];

(* Calculate the minimum distance for the trajectory of the nonlinear system that starts from initial point on the δ-diamond/simplex. *)
ini:=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==Apply[Sequence,{stabcond[pmin],pmin},2],List];firstit=1;

{dis,tdis,tsol}=compdisini[system,ini];gmindis=dis;

Print["NONLINEAR SYSTEM: ","    Minimum Distance from Initial Steady State = ",dis//N,"    Time of Closest Approach = ",tdis//N,"    Solution Time = ",tsol//N];

(* SOLVE - Solve for the global nonlinear saddle path *);
If[localmin==1,localmin=0;pmin=pmin+delta/100{1/3,1/3,-2/3}];

{gtmin,gminsol,pmin,epsilon}=compdis2[pmin,epsilon];

If[gtmin==(10^100),
Print["Either the system is highly ill-conditioned or you might have encountered a local minima. Perturing the center of the initial ring and resetting 'precision' to its initial value"];precision=precisionini;localmin=1;Goto[begin]];

(* ERROR CHECK - Check for the accuracy of the solution over the time period of interest. *) 
{firstsim,firstrun, delta, precision, doneflag}=errorCheck[];
If[doneflag==1,Goto[end]];

];
 
(* DURATION CHECK - Check the duration for which the solution is accurate within tolerance 'tol'.*);
Label[end];
durationCheck[];

{gtmin,gminsol}];

checkInput[___]:=Block[{},

If[Not[Length[sys]==Length[vars]==Length[varso]==Length[varsf]],
	Message[solve::sysLength,"'solve'","'sys'","'vars'","'varso'","'varsf'"];Abort[]];

If[Not[Length[vars]==Length[Cases[vars,_Symbol]]],Message[solve::vars,"'solve'","'vars'"];Abort[]];

If[Not[Length[varso]==Length[Cases[varso,_?NumericQ]]],Message[solve::origVals,"'solve'","'varso'"];Abort[]];

If[Not[Length[varsf]==Length[Cases[varsf,_?NumericQ]]],Message[solve::finVals,"'solve'","'varsf'"];Abort[]];

If[Not[Length[varsendo]==Length[Cases[varsendo,{_Symbol,_List}]]]||Not[Length[Flatten[varsendo]]==Length[Cases[Flatten[varsendo],_Symbol]]],
	Message[solve::endoVars,"'solve'","'varsendo'","'vars'"];Abort[];];

If[Not[Length[init]==Length[Cases[init,_?NumericQ]]],Message[solve::initVals,"'solve'","'init'"];Abort[]];

];

makeSys[___]:=Block[{system},

system=-sys/. Thread[MapThread[(#1)&,Transpose[varsendo]]->
					MapThread[ToExpression[StringJoin[ToString[ #1],"[",ToString[#2],"]"]] /.List->Sequence&, Transpose[varsendo]],
						List];

If[Length[Cases[system/.Thread[vars->varsf,List],_?NumberQ]]≠Length[vars],
	Message[solve::sysDESys,"'solve'","'sys'","'vars'","'varsendo'"];Abort[];];

system=Thread[Map[ToExpression[StringJoin[ToString[#],"'","[",ToString[t],"]"]]&,vars]==  
			(system=system/.Thread[vars->Map[ToExpression[StringJoin[ToString[#],"[",ToString[t],"]"]]&,vars],List]),
				List];

system];

linearize[___]:=Block[{ sysL,A,vals,vecs, negvalsPos, negvecsT},

dev= (10^(-50));

sysL= sys/.Thread[MapThread[(#1)&,Transpose[varsendo]]->
			MapThread[ToExpression[StringJoin[ToString[ #1],"[",ToString[#2],"]"]] /.List->Sequence&, Transpose[varsendo]],
				List];
A=Transpose[Table[(sysL/.Thread[vars->varsf+Table[If[j==i,dev,0],{j,1,Length[vars]}],List])-(sysL/.Thread[vars->varsf+Table[If[j==i,-dev,0],{j,1,Length[vars]}],List]),{i,1,Length[vars]}]/(2dev)];
{vals,vecs}=Eigensystem[A];
valratio = First[Take[Sort[vals],{2}]/Take[Sort[vals],{3}]];
Print[vals//N];Print[valratio//N];

negvalsPos=Flatten[Position[Map[Chop[Re[#]]<0&,vals],True]];
If[Length[negvalsPos]≠ 3,
	Message[solve::sysError,"'solve'"];
	Print["The eigenvalue of the linearized system are: "];
	Print[vals//N];Abort[]];
negvecsT = Transpose[Part[vecs,negvalsPos]];

sysLsol[x_]:=Evaluate[varsf-negvecsT.DiagonalMatrix[Map[Exp,Select[vals,Chop[Re[#]]<0&]x]].
						Inverse[Part[negvecsT,Table[i,{i,Length[vars]-2,Length[vars]}]]].Take[varsf-varso,-3]];

stabcond[{x_,y_,z_}]:=Evaluate[First[Simplify[Drop[vars,-3]/.
						Solve[(Drop[vars-varsf,-3])==Part[negvecsT,Table[i,{i,1,Length[vars]-3}]].
						Inverse[Part[negvecsT,Table[i,{i,Length[vars]-2,Length[vars]}]]].({x,y,z}-Take[varsf,-3]),
							Drop[vars,-3] ,WorkingPrecision->Infinity]]]];
A];

compdisini[system_,ini_]:=Block[{sol,tsol,dis,tdis},

sol=NDSolve[{system,ini},vars,{t,0,Tmax},
		Method->{"EventLocator", "Event"->box,"EventAction":> Throw[Null, "StopIntegration"],"EventLocationMethod"->"StepBegin"},
		MaxSteps->100000,WorkingPrecision->digits];

tsol=sol[[1,1,2,1,1,2]];{dis,tdis}=distance[sol,tsol];

{dis,tdis,tsol}];

distance[sol_,tsol_]:=Block[{dis,dislast,tdis,steps,tbeg,tend,intsz,traj,min},
If[(*firstit==1*)True,(* With True here, firstit has no effect on evaluation of minimum distance *)

	dis= (10^7);dislast=(10^8);steps=Max[4Floor[tsol],450] ;tbeg=0;tend=tsol;intsz=tend-tbeg;

	While[ (dislast-dis)>gmindis (10^(-7)),
		dislast=dis;
		traj=Table[Norm[Flatten[{vars[[-3]][t]-varso[[-3]],vars[[-2]][t]-varso[[-2]],vars[[-1]][t]-varso[[-1]]}/.sol],2],{t,tbeg,tend,(intsz/steps)}];
		dis=Min[traj];tdis=tbeg+(Position[traj,dis][[1,1]]-1)(intsz/steps);
		tbeg=Max[tdis-(intsz/steps),0];tend=Min[tdis+(intsz/steps),tsol];intsz=tend-tbeg;steps=40 ;
	],

	tbeg=0;tend=tsol;

	min=FindMinimum[First[Norm[Flatten[{vars[[-3]][t]-varso[[-3]],vars[[-2]][t]-varso[[-2]],vars[[-1]][t]-varso[[-1]]}]]/.sol],
					{t,9tsol/10,tbeg,tsol},WorkingPrecision->digits,AccuracyGoal->7-Log[10,gmindis]];
	{dis,tdis}={min[[1]],min[[2,1,2]]};];

{dis,tdis}];

(* Set up for solving for the global nonlinear saddle path *)
compdis2[p_,ϵ_]:=Block[{rmin,rtmin,rminsol,rpmin,rit,trajtrack},

{pmin,epsilon}={p,ϵ};rit=0;trajtrack=0;

{gmindis,gtmin,pmin,epsilon}=compdis1[pmin,epsilon];

trajtrack=0;rit=1;

While[True,
Label[incprec];
Print[" "];Print["BEGINNING SEARCH ON RING ",rit," WITH CENTER ",pmin];
{rmin,rtmin,rminsol,rpmin}=compdis[pmin,epsilon];

If[rtmin==(10^100)&& precision≤80,
	precision=precision+4;$MinPrecision=digits=precision;Print["Increased 'precision' by 4 digits to ",precision," digits."];Goto[incprec],
	If[rtmin==(10^100),gtmin=rtmin; Break[]]];
If[rmin≤gmindis,
	gmindis=rmin;gtmin=rtmin;gminsol=rminsol;pmin=rpmin;epsilon=epsilon/2,
	epsilon=epsilon/2];

Print["AT THE END OF SEARCH ON RING ",rit,":    CURRENT MINIMUM DISTANCE = ",gmindis//N,"     TIME OF CLOSEST APPROACH = ",gtmin];rit++;
If[gmindis<tol/3,
	Print["Target Tolerance = ",tol/3//N,"   Tolerance Achieved = ",gmindis]; Break[] ];
];

{gtmin,gminsol,pmin,epsilon}];

(* Set up for solving for the global nonlinear saddle path *)
compdis1[p_,ϵ_]:=Block[{irit,irmin,irtmin,irminsol,irpmin,srchmin,srchtmin,srchpmin},

irit=1;{pmin,epsilon}={p,ϵ};srchmin=gmindis;srchpmin=pmin;srchtmin=gtmin;

While[True,
{irmin,irtmin,irminsol,irpmin}=compdis[pmin,epsilon];Print["INITIAL RING ",irit,":    MINIMUM DISTANCE = ",irmin//N,"     TIME OF CLOSEST APPROACH = ",irtmin];irit++;
If[irmin≤srchmin,
	gmindis=srchmin=irmin;srchtmin=irtmin;srchpmin=irpmin;
	If[firstrun==1,epsilon=epsilon/2,epsilon=2epsilon(*/2*)],
	If[firstrun==1,epsilon=4epsilon,epsilon=(*4*)epsilon/2]; Break[] ];
];
Print["AFTER INITIAL SEARCH:    MINIMUM DISTANCE = ",srchmin//N,"     TIME OF CLOSEST APPROACH = ",srchtmin];

{srchmin,srchtmin, srchpmin,epsilon}];

(* Set up for solving for the global nonlinear saddle path *)
compdis[p_,ϵ_]:=Module[{solveit,theta,ntraj,pini,ini,sol,dis,tdis,j,jdis,jtdis,dissol,jp,tab,mindis,mindis2,j1,j2,ii,tmin,jmin,minsol,tmin2,pmin1,jmin2,minsol2,pmin2},

firstit=1;ntraj=ntrajini;

solveit:=(theta = 2π j/ntraj;pini=p+ϵ(Cos[theta]dir{1/(Sqrt[2]),-1/(Sqrt[2]),0}+ Sin[theta]dir{-1/(Sqrt[6]),-1/(Sqrt[6]),(Sqrt[2])/(Sqrt[3])});
ini=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==Apply[Sequence,{stabcond[pini],pini},2],List];
sol=NDSolve[{system,ini},vars,{t,0,Tmax},
		Method->{"EventLocator", "Event"->box,"EventAction":> Throw[Null, "StopIntegration"],"EventLocationMethod"->"StepBegin"},
			MaxSteps->Infinity,WorkingPrecision->digits,AccuracyGoal->Infinity];{dis,tdis}=distance[sol,sol[[1,1,2,1,1,2]]];);

For[j=1,j≤ ntraj  ,j++,
	If[trajtrack==0||(trajtrack==1 && (j==minpos[[1]]||j==minpos[[2]])),
		solveit;jdis[j]=dis;jtdis[j]=tdis;dissol[j]=sol;jp[j]=pini,jdis[j]=(10^7)];
];

tab=Table[jdis[Mod[ii-1,ntraj]+1]-jdis[Mod[ii,ntraj]+1],{ii,1,ntraj}];firstit=0;minpos=Flatten[Position[Table[tab[[Mod[ii-2,ntraj]+1]]>0 &&tab[[Mod[ii-1,ntraj]+1]]<0,{ii,1,ntraj}],True]];

If[rit≠0 &&Length[minpos]==2,trajtrack= 1];
If[rit≠ 0,
	Print["Local minima from initial search over the ring."];
	For[ii=1,
		ii≤Length[minpos],ii++,
		Print[ii,":" ,"   {mindis, tdis, jmin, ntraj} = ",{jdis[minpos[[ii]]],jtdis[minpos[[ii]]]//N,minpos[[ii]],ntraj}];
	];];

mindis=mindis2=(10^7);j1=1;j2= ntraj;

For[ii=1,ii≤Length[minpos],ii++,
	{mindis,tmin,jmin,ntraj,minsol}={jdis[minpos[[ii]]],jtdis[minpos[[ii]]],minpos[[ii]],ntrajini,dissol[minpos[[ii]]]};
	pmin1=jp[minpos[[ii]]];{j1,j2}=2{jmin-1,jmin+1};{ntraj,jmin}=2{ntraj,jmin};

	While[True,
	For[j=j1+1,j≤ j2  ,j+=2,
		solveit;
		If[ Max[Abs[Apply[(#1/#2-1)&,{pini,pmin1}]]]<(10^(-precision)),
			Print["You may be losing precision in calculations."] ;tmin2=(10^100);Goto[end2];];
		If[ Abs[mindis-dis]≤gmindis/100000,Goto[end1];];
		If[dis≤mindis,mindis=dis;tmin=tdis;jmin=j;minsol=sol;pmin1=pini;Break[];];
	];

	{j1,j2}=2{jmin-1,jmin+1};{ntraj,jmin}=2{ntraj,jmin};
	];

	Label[end1];

	If[mindis<mindis2,
	mindis2=mindis;tmin2=tmin;jmin2=jmin;minsol2=minsol;pmin2=pmin1];
	If[mindis2<gmindis,Break[]];

];

If[rit≠ 0 ,Print["Ring ",rit,":    Minimum Distance = ",mindis2//N,"    Time of Closest Approach = ",tmin2//N,"    No. of Subdivisions of ϵ-ring > ",jmin2];];

Label[end2];
Remove[solveit,theta,ntraj,pini,ini,sol,dis,tdis,j,jdis,jtdis,dissol,jp,tab,mindis,j1,j2,ii,tmin,jmin,minsol,pmin1,jmin2];
{mindis2,tmin2,minsol2,pmin2}];

errorCheck[___]:=Block[{diffprec,diffdel, satisprec, notsatisprec, notsatisdel, increaseprec,decreasedel, step2, step3, congratulations, done,gtmin1,gminsol1},

doneflag=0;

(* Messages *)
diffprec:=Print["Difference in values of variables due to latest increase in precision = ", sysDiff[0,T]//N];
diffdel:=Print["Difference in values of variables due to latest decrease in delta = ",sysDiff[0,T]//N];

satisprec= "THE PRECISION IS SATISFACTORY. The variables changed by less than 'tol' due to latest increase in 'precision' for current δ-diamond\\simplex.";
notsatisprec="The precision is NOT satisfactory. The variables changed by greater than 'tol' due to latest increase in 'precision' for current δ-diamond\\simplex.";
increaseprec:=Print["Increased 'precision' by 4 digits to ",precision," digits."];

notsatisdel="The δ-diamonds are NOT tight enough. The variables changed by greater than 'tol' due to latest decrease in 'delta'.";
decreasedel:=Print["Decreased 'delta' by a half to ",N[delta,precision],"."];

step2="Starting Step 2: Checking for the precision of the candidate solution.";step3="Starting Step 3: Checking for the accuracy of the candidate solution.";

congratulations= "CONGRATULATIONS! The variables changed by less than 'tol' due to latest decrease in 'delta'.";
done = "You are done! Hope to see you in another model. Continue with the next cell.";

(* Check candidate solution *)
Which[
firstsim == 1 && firstrun == 1, 
	{gminsol,gtmin,pmin,5epsilon}>>solution1;Print[step2];
	Print[satisprec];Print[step3];firstsim=0;firstrun=1;delta=2delta/3(*2*);decreasedel;(* This line leads to skipping of precision check. *)
	(*precision=precision+4;increaseprec;firstrun=0*),

firstsim == 1 && firstrun == 0,
	{gminsol1,gtmin1}=(<<solution1)[[{1,2}]];{gminsol,gtmin,pmin,5epsilon}>>solution1;diffprec;
	If[Max[sysDiff[0,T]]>tol,
		Print[notsatisprec];precision=precision+4;increaseprec,
		Print[satisprec];Print[step3];firstsim=0;firstrun=1;delta=2delta/3(*2*);decreasedel],

firstsim == 0  && firstrun == 1,
	{gminsol,gtmin,pmin,5epsilon}>>solution2;{gminsol1,gtmin1}=(<<solution1)[[{1,2}]];diffprec;
	If[Max[sysDiff[0,T]]>tol,
		{gminsol,gtmin,pmin,5epsilon}>>solution1;Print[notsatisdel];delta=2delta/3(*2*);decreasedel;(* This line leads to skipping of precision check. *)
		(*firstrun=0;precision=precision+4;increaseprec*),
		diffdel;Print[congratulations];
		Print[done];doneflag=1],

firstsim == 0  && firstrun == 0,
	{gminsol1,gtmin1}=(<<solution2)[[{1,2}]];{gminsol,gtmin,pmin,5epsilon}>>solution2;diffprec;
	If[Max[sysDiff[0,T]]>tol,
		Print[notsatisprec];precision=precision+4;increaseprec,
		Print[satisprec];{gminsol1,gtmin1}=(<<solution1)[[{1,2}]];diffdel;
		If[Max[sysDiff[0,T]]>tol,
			{gminsol,gtmin,pmin,5epsilon}>>solution1;Print[notsatisdel];firstrun=1;delta=delta/2;decreasedel,
			Print[congratulations];Print[done];doneflag=1];];
];   

{firstsim,firstrun,delta,precision,doneflag} ];

sysDiff[t1_?NumberQ,t2_?NumberQ]:=Table[Max[Abs[Table[((vars[[i]][gtmin-t]/.gminsol)-(vars[[i]][gtmin1-t]/.gminsol1))[[1]],{t,t1,t2,1/10}]]],{i,1,Length[vars]}];

durationCheck[___]:=Block[{gminsol1,gtmin1,j},

{gminsol1,gtmin1}=(<<solution1)[[{1,2}]]; {gminsol,gtmin}=(<<solution2)[[{1,2}]];

For[j=T+1,j≤IntegerPart[gtmin1],j++,
	If[Max[sysDiff[j-1,j]]>tol,
		Print["The solution is accurate within tolerance 'tol' = ",tol//N," upto ",j-1," years."];Break[]];];

If[j==IntegerPart[gtmin1]+1,Print["The solution is accurate within tolerance 'tol' = ",tol//N," upto ",j-1," years."];];
Print["The solution upto ",T," years is accurate to within ",Max[sysDiff[0,T]]//N,". Continue with next cell."];

];

End[];

(* Epilog *)
Protect[Evaluate[$Context<>"*"]];